'''
nodalDiff.py - Paul Talbot, NE553, Oregon State University, Spring 2011
Used after runProblem.py to solve diffusion-based calculations.  Solves 
matrix operations of the form Mx=C, where M is a 10x10 matrix, x is a 
vector of scalar flux values, and C is a solution vector.  Uses Gauss-
Seidel methods to invert and solve for x.
'''

import gaussSeidel as GS
import time
import buildMesh as BM
import runProblem as RP


case1='U1'
case2='M1'
abls=[BM.Mesh(case1,'s10'),BM.Mesh(case2,'s10')]
for ably in abls:
  RP.runProblem(mesh=ably)

converged=False
  #guess F, k values
while not converged:
  #do 10x10 fast


def buildM(mesh,g): #10x10 coefficient matrix
  a=mesh.ablys[0].groups[g-1] #shorthand
  b=mesh.ablys[1].groups[g-1] #same
  adx=mesh.ablys[0].width
  bdx=mesh.ablys[1].width
  M=[
      [0,1,-3,6,-10,  0,0,0,0,0], #1
      [0,0,0,0,0,     0,1,3,6,10], #2
      [0,a.D/adx,3*a.D/(2*adx),3*a.D/(8*adx),-25*a.D/(64*adx),
          0,-b.D/bdx,-3*b.D/(2*bdx),-3*b.D/(8*bdx),25*b.D/(64*bdx)], #3
      [a.fR,a.fR,a.fR,a.fR,a.fR,  b.fL,b.fL,b.fL,b.fL,b.fL], #4
      [adx*a.sigr,0,-12*a.D/adx,0,-40*a.D/adx,  0,0,0,0,0], #5
      [0,0,0,0,0,  bdx*b.sigr,0,-12*a.D/bdx,0,-40*b.D/bdx], #6
      [0,a.sigr*adx/3.,0,-20*a.D/adx,0,  0,0,0,0,0], #7
      [0,0,0,0,0,  0,b.sigr*bdx/3,0,-20*b.D/bdx,0], #8
      [0,0,adx*a.sigr/5.,0,-28*a.D/adx,  0,0,0,0,0], #9
      [0,0,0,0,0,  0,0,bdx*b.sigr/5.,0,-28*b.D/bdx] #10
    ]
  #for i,row in enumerate(M):
  #  print 'Row:',i,' Size:',len(row)
  return M

def buildC(mesh,g):
  a=mesh.ablys[0].groups[g-1] #shorthand
  b=mesh.ablys[1].groups[g-1]
  adx=mesh.ablys[0].width
  bdx=mesh.ablys[1].width
  if g==1: #fast case
    Sa=Sb=mesh.Fco
    ka=kb=1/mesh.K
  else: #thermal case
    Sa=a.phiCo #vector of Legendre Coefficients
    Sb=b.phiCo #same
    ka=a.sigr #removal coefficient
    kb=b.sigr
  C=[[
    0, #1
    0, #2
    0, #3
    0, #4
    adx*ka*Sa[0], #5
    bdx*kb*Sb[0], #6
    -adx*ka*Sa[1]/3., #7
    -bdx*kb*Sb[1]/3., #8
    adx*ka*Sa[2]/5., #9
    bdx*kb*Sb[2]/5. #10
    ]]
  return C

def runDiffusion(mesh):
  init=time.clock()
  M1=buildM(mesh,1)
  M2=buildM(mesh,2)
  C1=buildC(mesh,1)
  C2=buildC(mesh,2)
  elapsed=time.clock()-init
  phi1=GS.doMatrixSolve(M1,C1)
  phi2=GS.doMatrixSolve(M2,C2)
  print phi1
  print phi2
  print 'Diffusion run time:',elapsed


